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ABSTRACT 

We present results of general relativistic simulations of collapsing supermassive stars with and 
without rotation using the two-dimensional general relativistic numerical code Nada, which solves 
the Einstein equations written in the BSSN formalism and the general relativistic hydrodynamics 
equations with high resolution shock capturing schemes. These numerical simulations use an equation 
of state which includes effects of gas pressure, and in a tabulated form those associated with radiation 
and the electron-positron pairs. We also take into account the effect of thermonuclear energy released 
by hydrogen and helium burning. We find that objects with a mass of m 5 x 10 5 Mq and an initial 
metallicity greater than Zcno ~ 0.007 do explode if non-rotating, while the threshold metallicity for 
an explosion is reduced to Zcno ~ 0.001 for objects uniformly rotating. The critical initial metallicity 
for a thermonuclear explosion increases for stars with mass w 10 6 M Q . For those stars that do not 
explode we follow the evolution beyond the phase of black hole formation. We compute the neutrino 
energy loss rates due to several processes that may be relevant during the gravitational collapse of 
these objects. The peak luminosities of neutrinos and antineutrinos of all flavors for models collapsing 
to a BH are L„ ~ 10 55 erg/s. The total radiated energy in neutrinos varies between E u ~ 10 56 ergs 
for models collapsing to a BH, and E v ~ 10 45 — 10 46 ergs for models exploding. 
Subject headings: Supermassive stars 



1. INTRODUCTION 

There is large observational evidence of the presence 
of supermassive black holes (SM BH) in the centres of 
most nearby galaxies (I Reed (19981 ) . The dynamical evi- 
dence related to the orbital motion of stars in the cluster 
surrounding Sgr A* indic ates the presence o f a SMBH 
with mass « 4 x 1O 6 M (|Genzel et all 120001 ). In addi- 
tion, the observed correlation between the central black 
hole masses and the stellar velocity dispersion of the 
bulge of the host galaxies suggests a direct connection 
between the formation and evolution of galaxies and 
SMBH (IKormendv fe Gebhardt ll200l . 

The observation of luminous quasars detected at red- 
shifts higher than 6 in the Sloan Digital Sky Survey 
(SDSS) implies that SMBH with masses ~ 10 9 M o , which 
are believed to be the engines of such powerful quasars, 
were formed within the first billion years after the Big 
Bang (e.g. Fan 2006 for a recent review). However, it is 
still an open question how SMBH seeds form and grow 
to re ach such h igh masses in such a short amount of 
time (lReesll2001l ). 

A number of different routes based on stellar dynamical 
processes, hydrodynamical processes or a combination of 
both have been suggested (e.g. Volonteri 2010 for a re- 
cent review) . One of the theoretical scenarios for SMBH 
seed formation is the gravitational collapse of the first 
generation of stars (Population III stars) with masses 
M ~ IOOMq that are expected to form in halos with 
virial temperature T V i r < lO^K at z ~ 20 — 50 where 
cooling by molecular hydrogen is effective. As a result 
of the gravitational collapse of such Pop III stars, very 
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massive B Hs would form and then grow via merger and 
accretion (lHaiman fc L oeb 2001; lYoo fc Miralda-Escudl 
l2lMlATvarez et al.M2009lL 

Another possible scenario proposes that if sufficient 
primordial gas in massive halos, with mass ~ 10 8 M Q , is 
unable to cool below TU r > 10 4 7^, it may lead to the for- 
mation of a supermas sive object (|Bromm fc Loebl feOOS: 
iBegelman et all 120061) . which would eventually collapse 
to form a SMBH. This route assumes that fragmenta- 
tion, which depends on efficient cooling, is suppressed, 
possibly by the presence of sufficiently strong UV radi- 
ation, that prevents the formation of molecular hydro- 
gen in an environ ment with metallicity smaller than a 
given critical value (jSantoro fc S hull 2006S; lOmukai et al.l 
2008). Furthermore, fragmentation may depend on 
the turbulence present within the inflow of gas, and 
on t he mechanism redistributing its angular momen- 
tum (jBcgelman & Shlo smanl 12009 ). T he "bars-within- 

bars" mechanism (fShlosma n et al.lll9 89: Beac lman et al.l 
2006) is a self-regulating route to redistribute angular 
momentum and sustain turbulence such that the inflow 
of gas can proceed without fragmenting as it collapses 
even in a metal-enriched environment. 

Depending on the rate and efficiency of the inflow- 
ing mass, there may be different outcomes. A low 
rate of mass accumulation would favor the formation 
of isentropic supermassive stars (SMS), with mass > 
5 x 10 4 M Q , which then would evolve as equili brium con- 
figurations dominated by radiation pressure (Ibcn 1963; 
iHovle fc Fowler|[l963t lFowler|[l964|) . A different outcome 
could result if the accumulation of gas is fast enough so 
that the outer layers of SMS are not thermally relaxed 
during much of their lifetime , thus having an entropy 
stratification (Bcgclman 2009). 

A more exotic mechanism that could eventually lead to 
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a SMS collapsing into a SMBH is the formation and evo- 
lution of supermassive dark matter stars (SDMS) (Spol- 
yar et al 2008). Such stars would be composed pri- 
marily of hydrogen and helium with only about 0.1% 
of their mass in the form of dark matter, however they 
would shine due to dark matter annihilation. It has re- 
cently bee n pointed out that SDMSs could reach masses 
- 10 5 M Q (jFreese et al.ll20Tol) . Once SDMSs run out of 
their dark matter supply, they experience a contraction 
phase that increases their baryon density and tempera- 
ture, leading to an environment where nuclear burning 
may become important for the subsequent stellar evolu- 
tion. 

If isentropic SMS form, their quasi-stationary evolu- 
tion of cooling and contraction will drive the stars to 
the onset of a general relativistic gravit ational instabil- 
ity le a ding to their gr avitational collapse ( Chandra sekharl 
119641: iFowlerl 119641) . and possibly also to the for- 
mation of a SMBH. The first numerical simula- 
tions, within the p o st-Ne wtonian approximation, of 
lAppenzeller fc Frickel ([1972) concluded that for spheri- 
cal stars with masses greater than 1O 6 M thermonuclear 
reactions have no major effect on the collapse, while less 
massive stars exploded due to hydrogen burning. Later 
IShapiro fc Teukolskvl ()1979l ) performed the first relativis- 
tic simulations of the collapse of a SMS in spherical 
symmetry. They were able to follow the evolution un- 
til the formation of a BH, although their investi g ations 
did not include any microphysics. iFuller et all (1986) 
revisited the work of lAppenzeller fc Frickel (|1972 ) and 
performed simulations of non-rotating SMS in the range 
of 10 5 -10 6 M Q with post-Newtonian corrections and de- 
tailed microphysics that took into account an equation of 
state (EOS) including electron-positron pairs, and a re- 
action network describing hydrogen burning by the CNO 
cycle and its break-out via the rapid proton capture (rp)- 
process. They found that SMS with zero initial metallic- 
ity do not explode, while SMS with masses larger than 
10 5 M© and with metallicity Zn no > 0.005 do explode. 

More recently iLinke et al.l (|2001l ) carried out general 
relativistic hydrodynamic simulations of the spherically 
symmetric gravitational collapse of SMS adopting a 
spacetime foliation with outgoing null hypersurfaces to 
solve the system of Einstein and fluid equations. They 
performed simulations of spherical SMS with masses in 
the range of 5 x 10 5 M Q - 1O 9 M using an EOS that 
accounts for contributions from baryonic gases, and in 
a tabulated form, radiation and electron-positron pairs. 
They were able to follow the collapse from the onset 
of the instability until the point of BH formation, and 
showed that an apparent horizon (AH) enclosing about 
25% of the stellar material was formed in all cases when 
si mulations stopped. 

Shibata & Shaping (|2002l ) carried out general relativis- 
tic numerical simulations in axisymmetry of the collapse 
of uniformly rotating SMS to BHs. They did not take 
into account thermonuclear burning, and adopted a T- 
law EOS, P = (T - l)pe with adiabatic index T = 4/3, 
where P is the pressure, p the rest-mass density, and 
e the specific internal energy. Although their simula- 
tions stopped before the final equilibrium was reached, 
the BH growth was followed until about 60% of the mass 
had been swallowed by the SMBH. They estimated that 



about 90% of the total mass would end up in the final 
SMBH with a spin parameter of J/M 2 ~ 0.75. 

The gravitational collapse of differentially rotat- 
ing SMS in t hree dimensions was investigated by 
Saiio & Hawkej <|2009f ) . who focused on the post-BH evo- 
lution, and also on the gravitational wave (GW) signal re- 
sulting from the newly formed SMBH and the surround- 
ing disk. The GW signal is expected to be emitted in 
the low frequency LISA band (10~ 4 - lO" 1 Hz). 

Despite the progress made, the final fate rotating isen- 
tropic SMS is still unclear. In particular, it is still an 
open question for which initial mctallicities hydrogen 
burning by the /3-limited hot CNO cycle and its break- 
out via the 15 0(a, 7) 19 Ne reaction (rp-process) can halt 
the gravitational collapse of rotating SMS and generate 
enough thermal energy to lead to an explosion. To ad- 
dress this issue, we perform a series of general relativis- 
tic hydrodynamic simulations with a microphysical EOS 
accounting for contributions from radiation, electron- 
positron pairs, and baryonic matter, and taking into ac- 
count the net thermonuclear energy released by the nu- 
clear reactions involved in hydrogen burning through the 
pp-chain, cold and hot CNO cycles and their break-out 
by the rp-process, and helium burning through the 3- 
a reaction. The nu merical simulations we re carried out 
with the Nada code (|Montero et al J 120081) . which solves 
the Einstein equations coupled to the general relativistic 
hydrodynamics equations. 

Greek indices run from to 3, Latin indices from 1 to 3, 
and we adopt the standard convention for the summation 
over repeated indices. Unless otherwise stated we use 
units in which c = G = 1 . 

2. BASIC EQUATIONS 

Next we briefly describe how the system of Einstein 
and hydrodynamic eq uations are implemen ted in the 
Nada code. We refer to iMontero et all (|2008l ) for a more 
detailed description of the main equations and thorough 
testing of the code (namely single BH evolutions, shock 
tubes, evolutions of both spherical and rotating relativis- 
tic stars, gravitational collapse to a BH of a marginally 
stable spherical star, and simulations of a system formed 
by a BH surrounded by a self-gravitating torus in equi- 
librium). 

2.1. Formulation of Einstein equations 
2.1.1. BSSN formulation 

We follow the 3+1 formulation in which the spacetime 
is foliated into a set of non-intersecting spacelike hyper- 
surfaces. In this approach, the line element is written in 
the following form 



ds 2 



-(a 2 - PiP^dt 2 + 2/3 i dx i dt + ^dx l dx\ (1) 



where a, f3 l and jij are the lapse function, the shift three- 
vector, and the three-metric, respectively. The latter is 
defined by 

lnv= g^+n^riu, (2) 

where is a timelike unit-normal vector orthogonal to 
a spacelike hypersurface. 

We make use of the BSSN formulation (|Nakamura 
1987; Shibata & Nakamura 1995; Baumgartc & Shapiro 
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1999) to solve the Einstein equations. Initially, a confor- 
mal factor ef> is introduced, and the conformally related 
metric is written as 



(3) 



such that the determinant of the conformal metric, 7^-, 



is unity and 



ln(7)/12, where 7 = dct^. 



We 



also define the conformally related traceless part of the 
extrinsic curvature Ka, 



1 



•7ij 



K 



(4) 



where K is the trace of the extrinsic curvature. We 
evolve the conformal factor defined as x = e -4 "^ 
(|Campanelli et aill2006() . and the auxiliary variables T l , 
known as the conformal connection functions, defined as 



-jkfi 
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(5) 



where T r . fe are the connection coefficients associated 
with jij. 

During the evolution we also enforce the constraints 
Tr(Aij) = and det(7y ) = 1 a t every time step . 

We use the Cartoon method (jAlcubierre et al.|[200lD to 
impose axisymmetry while using Cartesian coordinates. 

2.1.2. Gauge choices 

In addition to the BSSN spacetime variables 
(jij, Aij, K, x, T 1 ), there are two more quantities left un- 
determined, the lapse, a and the shift vector j3 l . We used 
the s o-called "non-advective 1+log" slicing (|Bona et al.l 
1997), by dropping the advective term in the "1+log" 
slicing condition. In this case, the slicing condition takes 
the form 

d t a = -2aK. (6) 

For the shift vector, we choose t he "Gamma-freezing con- 
dition" (jAlcubierre et alJl2003t > written as 



d t p = -b\ 
d t B i = d t r -nB\ 



(7) 
(8) 



where n is a constant that acts as a damping term, orig- 
inally introduced both to prevent long term drift of the 
metric functions and to prevent oscillations of the shift 
vector. 

2.2. Formulation of the hydrodynamics equations 

The general relativistic hydrodynamics equations, ex- 
pressed through the conservation equations for the stress- 
energy tensor T^ v and the continuity equation are 



V /i T' u ' 







V M (pu») = 0, 



(9) 



where p is the rest-mass density, u M is the fluid four- 
velocity and V is the covariant derivative with respect to 
the spacetime metric. Following Shibatal •: 2()(K|). the gen- 
eral relativistic hydrodynamics equations are written in 
a conservative form in cylindrical coordinates. Since the 
Einstein equations are solved only in the y = plane with 
Cartesian coordinates (2D), the hydrodynamic equations 



are rewritten in Cartesian coordinates for y = 0. The fol- 
lowing definitions for the hydrodynamical variables are 
used 

p, = pWe 6 *, (10) 



ir hW 
u 4 = hui, 

e 6<t> p 

e= — T^n" = hW - — , 
P* pW 

W = au\ 



(11) 
(12) 

(13) 
(14) 

where W and h are the Lorentz factor and the specific 
fluid enthalpy respectively, and P is the pressure. The 
conserv ed variables ar e p*, Ji = p*iii, £* = p*e. We 
refer to IShibatal (|2003f ) for further details. 

3. SUPERMASSIVE STARS AND MICROPHYSICS 

3.1. Properties of SMS 

Isentropic SMS are self-gravitating equilibrium config- 
urations of masses in the range of 10 4 — 10 8 M Q , which 
are mainly supported by radiation pressure, while the 
pressure of electron-positron pairs and of the baryon gas 
are only minor contributions to the EOS. Such configu- 
rations are well described by Newtonian polytropes with 
polytropic index n — 3 (adiabatic index T = 4/3). The 
ratio of gas pressure to t he total pressure (0) fo r spherical 
SMS can be written as (|Fowler fc Hovlelll966ft 



3l 

Ptot 



4.3 

A* 



Mo 
M 



1/2 



(15) 



where p is the mean molecular weight. Thus j3 « 10 2 
for M « 1O 6 M . 

Since nuclear burning timescales are too long for 
M > 1O 4 M , evolution of SMS proceeds on the Kelvin- 
Helmholtz timescale and is driven by the loss of energy 
and entropy by radiation as well as loss of angular mo- 
mentum via mass shedding in the case of rotating con- 
figurations. 

Although corrections due to the nonrelativistic gas of 
baryons and electrons and general relativistic effects are 
small, they cannot be neglected for the evolution. Firstly, 
gas corrections raise the adiabatic index slightly above 
4/3 



r * t + f + ° (/?2) 



(16) 



Secondly, general relativistic corrections lead to the 
existence of a maximum for the equilibrium mass as a 
function of the central density. For spherical SMS this 
means that for a given mass the star evolves to a critical 
density beyond whic h it is dynamically un stable against 
radial perturbations (Chandrasckhar 1964): 



Pcrit = 1.994 x 10 1 



0.5 



Mq 

M 



7/2 



gem 



-3 



(17) 



The onset of the instability also corresponds to a criti- 
cal value of the adiabatic index T cr n, i.e. configurations 



4 



Montero, Janka and Miiller 



become unstable when the adiabatic index drops below 
the critical value 



_ 4 2GM 

Trr-it — ~T~ 1.12 — . 

cnt 3 Rc 2 



(18) 



This happens when the stabilizing gas contribution to 
the EOS does not raise the adiabatic index above 4/3 to 
compensate for the destabilizing effect of general relativ- 
ity expressed by the second term on the righ-hand-side 
ofEq. (pj). 

Rotation can stabilize configurations against the ra- 
dial instability. The stability of rotating SMS with 
uniform rotation was analyzed by iBaumgarte &: Shapiro! 
(1999a. b). They found that stars at the onset of the 
instability have an equatorial radius R ~ 640GAf/c 2 , a 
spin parameter q = cJ/GM 2 « 0.97, and a ratio of rota- 
tional kinetic energy to the gravitational binding energy 
of T/Wk 0.009. 

3.2. Equation of State 

To close the system of hydrodynamic equations (Eq. [5]) 
we need to define the EOS. We follow a treatment which 
includes separately the baryon contribution on the one 
hand, and photons and electron-positron pairs contri- 
butions, in a tabulated form, on the other hand. The 
baryon contribution is given by the analytic expressions 
for the pressure and specific internal energy 



Ph = 



tb = 



KpT 

Mb 

2 KpT 
Pb 



(19) 



(20) 



where K the universal gas constant, T the temperature, 
e;, the baryon specific internal energy, and p^ is the mean 
molecular weight due to ions, which can be expressed as 
a function of the mass fractions of hydrogen (A), helium 
(Y) and heavier elements (metals) (Zcno) as 



1 

fib 



X 



Y Zcno 
4 + T4T 



(21) 



where (A) is the average atomic mass of the heavy ele- 
ments. We assume that the composition of SMS (approx- 
imately that of primordial gas) has a mass fraction of 
hydrogen X = 0.75 — Zcno and helium Y — 0.25, where 
the metallity Zcno — 1 — X — Y is an initial parameter, 
typically of the order of Zcno ~ 10~ 3 (see Table 1 de- 
tails). Thus, for the initial compositions that we consider 
the mean molecular weight of baryons is pi, rs 1.23 (i.e. 
corresponding to a molecular weight for both ions and 
electrons of p re 0.59). 

Effects associated with photons and the creation of 
electron-positron pairs are taken into account employ- 
ing a tabulated EOS. At temperatures above 10 9 K, not 
all the energy is used to increase the temperature and 
pressure, but part of the photon energy is used to create 
the rest-mass of the electron-positron pairs. As a result 
of pair creation, the adiabatic index of the star decreases, 
which means that the stability of the star is reduced. 

Given the specific internal energy, e and rest-mass den- 
sity, p, as evolved by the hydrodynamic equations, it is 



possible to compute the temperature T by a Newton- 
Raphson algorithm that solves the equation e* (p, T) = e 
for T, 



T n+ i = T n - (e*(p,T n ) - e) 



where n is the iteration counter. 



de*{p,T) 
df 



(22) 



3.3. Nuclear burning 

In order to avoid the small time steps, and CPU-time 
demands connected with the solution of a nuclear re- 
action network coupled to the hydrodynamic evolution, 
we apply an approximate method to take into account 
the basic effects of nuclear burning on the dynamics of 
the collapsing SMS. We compute the nuclear energy re- 
lease rates by hydrogen burning (through the pp-chain, 
cold and hot CNO cycles, and their break-out by the rp- 
process) and helium burning (through the 3-a reaction) 
as a function of rest-mass density, temperature and mass 
fractions of hydrogen X, helium Y and CNO metallicity 
Zcno- These nuclear energy generation rates are added 
as a source term on the right-hand-side of the evolution 
equation for the conserved quantity E*. 

The change rates of the energy density due to nuclear 
reactions, in the fluid frame, expressed in units of [erg 
cm~ 3 s _1 ] are given by: 

• pp-chain (|Clavtonlll98l : 



de 

— ) = p(2.38 x 10Vffii^T 6 
' pp 

e -33.80/T 6 °- 3333 -j 



where T 6 = T/10 6 K, and gn is given by 

g n = l + 0.0123T 6 03333 + 0.0109T 6 - 66666 - 
0.0009T 6 . 

3-a (fWiescher et al.lll999l) : 

=p(5.1xl0V^ 9 - 3 e- 4 ^), 
where T 9 = T/10 9 K. 

Cold-CNO cycle (|Shen fc Bildstenll2007l) : 



(23) 



(24) 



(25) 



= 4.4 x 10 25 p 2 XZ C no 



ccno 



(r -2/3 e -l5.231/T 9 1 / 3 + 



8.3 x 10-°T< 



5 T -3/2 -3.0057/Tg 



Hot-CNO cycle (IWiescher et aljfl999h : 



4.6 x l0 15 pZ CN o 



(26) 



(27) 



HCNO 
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• rp-process (jWiescher et al.lll999f) : 




= p(1.77 x 10 ie pYZ CN o 

rp 



29.96T 9 3 / 2 e -5-85/T 9 ^ ^ 

Since we follow a single fluid approach, in which we 
solve only the hydrodynamics equations Eq. (|9]) (i.e. we 
do not solve additional advection equations for the abun- 
dances of hydrogen, helium and metals), the elemental 
abundances during the time evolution are fixed. Nev- 
ertheless, this assumption most possibly does not af- 
fect significantly the estimate of the threshold metallicity 
needed to produce a thermal bounce in collapsing SMS. 
The average energy release through the 3-a reaction is 
about 7.275 MeV for each 12 C nucleus formed. Since the 
total energy due to helium burning for exploding mod- 
els is ~ 10 45 ergs (e.g. 9.0 x 10 44 ergs for model Six); 
and even considering that this energy is released mostly 
in a centr al region of the SM S containing 1O 4 M of its 
rest-mass (Fuller et al. 1986), it is easy to show that the 
change in the metallicity is of the order of 10~ n . There- 
fore, the increase of the metallicity in models experienc- 
ing a thermal bounce is much smaller than the critical 
metallicities needed to trigger the explosions. Similarly, 
the average change in the mass fraction of hydrogen due 
to the cold and hot CNO cycles is expected to be ~ 10% 
for exploding models. 

3.4. Recovery of the primitive variables 

After each time iteration the conserved variables 
(i.e. p*, J x , J y , J z , E*) are updated and the primitive hy- 
drodynamical variables (i.e. p, v x ,v v ,v z , e) have to be re- 
covered. The recovery is done in such a way that it allows 
for the use of a general EOS of the form P = P(p,e). 
We calculate a function f(P*) = P(p*,e*) - P* , where 
p* and e* depend only on the conserved quantities and 
the pressure guess P*. The new pressure is computed 
then iteratively by a Newton-Raphson method until the 
desired convergence is achieved. 

3.5. Energy loss by neutrino emission 

The EOS allows us to compute the neutrino losses due 
to the following processes, which become most relevant 
just before BH formation: 

• Pair annihilation (e + + e~ — > v + v): most impor- 
tant process above 10 9 K. Due to the large mean 
free path of neutrinos in the stellar medium at the 
densities of SMS the energy loss by neutrinos can 
be significant. For a 10 6 M Q SMS most of the en- 
ergy release in the form of neutrinos originates from 
this process. The rate s are computed u sing the fit- 
ting formula given bv lltoh et ail (|1996[ ). 

• Photo-neutrino emission (7 + e — > e^ + D + v): 
dominates at low tempe ratures T <■ 4 x 10 8 K and 
densities p < 10 5 gcm -3 (jltoh et alJll996D . 

• Plasmon decay (7 — >■ v + v): This is the least rele- 
vant process for the conditions encountered by the 
models we have considered because its importance 
increases at higher densities than those present in 



SMS. The rat es are computed u sing the fitting for- 
mula given bv lHaft et at] (| 19941) . 

4. COMPUTATIONAL SETUP 

The evolution equations are integrated by the method 
of lines, for which we use an optimal strongly stability- 
preserving (SS P) Runge-Kutta algori thm of fourth-order 
with 5 stages (jSpiteri fe Ruuthll2002l ). We use a second- 
order slope limiter reconstruction scheme (MC limitcr) to 
obtain the left and right states of the primitive variables 
at eac h cell interface, and a HLLE appro ximate Riemann 
solver (jHarten et al.lll983HEinfeldtlll988l) to compute the 
numerical fluxes in the x and z directions. 

Derivative terms in the spacetime evolution equa- 
tions are represented by a fourth-order centered finite- 
difference approximation on a uniform Cartesian grid ex- 
cept for the advection terms (terms formally like fi l diu), 
for which an upwind scheme is used. 

The computational domain is defined as < x < L and 
< z < L, where L refers to the location of the outer 
boundaries. We used a cell-centered Cartesian grid to 
avoid that the location of the BH singularity coincides 
with a grid point. 

4.1. Regridding 

Since it is not possible to follow the gravitational 
collapse of a SMS from the early stages to the phase 
of black hole formation with a uniform Cartesian 
grid (the necessary fine zoning would be computation- 
ally too demanding), we adopt a regridding procedure 
(jShibata fe Sha piro 2002). During the initial phase of 
the collapse we rezone the computational domain by 
moving the outer boundary inward, decreasing the grid 
spacing while keeping the initial number of grid points 
fixed. Initially we use NxN = 400 x 400 grid points, and 
place the outer boundary at L « 1.5r c where r e is the 
equatorial radius of the star. Rezoning onto the new grid 
is done using a polynomial interpolation. We repeat this 
procedure 3-4 times until the collapse timescale in the 
central region is much shorter than in the outer parts. 
At this point, we both decrease the grid spacing and also 
increase the number of grid points N in dependence of 
the lapse function typically as follows: NxN — 800 x 800 
if 0.8 > a > 0.6, N x N = 1200 x 1200 if 0.6 > a > 0.4, 
and N x N = 1800 x 1800 if a < 0.4. This procedure en- 
sures the error in the conservation of the total rest-mass 
to be less than 2% on the finest computational domain. 

4.2. Hydro- Excision 

To deal with the spacetime singularity from the newly 
formed BH we use the method of excising the matter 
content in a region within the horizon as proposed by 
IHawke et al.l ([2005) once an AH is found. This excision 
is done only for the hydrodynamical variables, and the 
coordinate radius of the excised region is allowed to in- 
crease in time. On the other hand, we do neither use 
excision nor artificial dissipation terms for the spacetime 
evolution, and solely rely on the gauge conditions. 

4.3. Definitions 

Here we define some of the quantities listed in Table 1. 
We compute the total rest-mass M* and the ADM mass 
M as 

pL p L 

M* = 47r / xdx \ p*dz, (29) 
Jo Jo 
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TABLE 1 

Main properties of the initial models studied. From left to right the columns show: model, gravitational mass, initial 
central rest-mass density, tfc/|w|, angular velocity on the equatorial plane at the surface, initial central temperature, 
metallicity, the fate of the star, radial kinetic energy after thermal bounce, and total neutrino energy output. 



Model 


M 


Pc 


T k /\W\ 


n 


T c 


Initial mctallicity 


Fate 




E v 




[1O 5 M ] 


[10- 2 g/cm 3 ] 




[10" 5 rad/s] 


[10 7 K] 


[10- 3 ] 




[10 56 erg] 


[erg] 


SI. a 


5 


2.4 








5.8 


5 


BH 




3.4 x 10 56 


Sl.b 


5 


2.4 








5.8 


6 


BH 






Sl.c 


5 


2.4 








5.8 


7 


Explosion 


5.5 


9.4 x 10 45 


R1.0 


5 


40 


0.0088 


2.49 


13 





BH 






Rl.a 


5 


40 


0.0088 


2.49 


13 


0.5 


BH 




5.4 x 10 56 


Rl.b 


5 


10 


0.0088 


2.49 


13 


0.8 


BH 






Rl.c 


5 


10 


0.0088 


2.49 


13 


1 


Explosion 


1.0 




Rl.d 


5 


40 


0.0088 


2.49 


13 


2 


Explosion 


1.9 


8.9 x 10 45 


S2.a 


10 


0.23 








2.6 


30 


BH 




6.8 x 10 56 


S2.b 


10 


0.23 








2.6 


50 


Explosion 


35 


8.0 x 10 46 


R2.a 


10 


12 


0.0087 


1.47 


9.7 


0.5 


BH 




3.1 x 10 56 


R2.b 


10 


12 


0.0087 


1.47 


9.7 


0.8 


BH 






R2.c 


10 


12 


0.0087 


1.47 


9.7 


1.0 


BH 






R2.d 


10 


12 


0.0087 


1.47 


9.7 


1.5 


Explosion 


1.5 


2.1 x 10 46 


Dl 


5 


6.9 x 10 4 


0.089 


540 


140 





BH 






D2 


6 


1.3 x 10 5 


0.128 


700 


170 





Stable/BH^, 







a If the contribution of pairs is not taken into account in the EOS (e.g., in the case T-law EOS or an EOS that includes only the 
radiation and pressure contributions in an analytic form) the model is stable against gravitational collapse. However, if the effect of e^ 
pairs is considered, the star becomes unstable against gravitational collapse due to the reduction of the adiabatic index associated with the 
pair creation. 



M = -2 



L rL 

xdx I dz 
ii Jo 



-2nEe 



AijA 



50 



-R 



-K' 



(30) 



where E = n ll n v T^ v (n M being the unit normal to the 

hypersurface) and R is the scalar curvature associated to 
the conformal metric 7^ . 

The rotational kinetic energy Tk and the gravitational 
potential energy W are given by 

pL p L 

Tk — 2n / x 2 dx I p r u y fldz, (31) 
Jo Jo 

where Q is the angular velocity. 

W = M - (M* + T k + E mt ) , (32) 
where the internal energy is computed as 

pL p L 

Eint ~4:n xdx / p*edz. (33) 
Jo Jo 

In axisymmetry the AH equation becomes a nonlinear 
ordinary d ifferential equation for the AH sh ape function, 
h = h{0) (|Shibatalll997t IThornburd[2007l) . We employ 
an AH finder that solves this ODE by a shooting method 
using dgh{6 = 0) = and dgh(8 = tt/2) = as boundary 
conditions. We define the mass of the AH as 



A 



(34) 



Mah = a/— , 

where A is the area of the AH. 

5. INITIAL MODELS 

The initial SMS are set up as isentropic objects. All 
models, except model D2, are chosen such that they are 
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T-Law EOS (r=4/3) 
EOS without e" and 
EOS with e* pairs 




FIG. 1. 



t[10 5 s] 



Time evolution of the central rest-mass density for 



i5tv 



model R1.0 (a uniformly rotating star with a mass M = 5 X 10 Mq 
with zero metallicity) for three different EOS (r-law and the mi- 
crophysical EOS with and without the electron-positron pair cre- 
ation). 



gravitationally unstable, and therefore their central rest- 
mass density is slightly larger than the critical central 
density required for the onset of the collapse of a configu- 
ration with given mass and entropy. A list of the different 
SMS we have considered is provided in Table 1. Models 
SI and S2 represent a spherically symmetric, nonrotat- 
ing SMS with gravitational mass of M = 5 x 10 5 M Q 
and M = 1 x 10 6 M Q , respectively, while models Rl 
and R2 are uniformly rotating initial models again with 
masses of M = 5 x 1O 5 M and M = 1 X 10 6 M Q , re- 
spectively. The rigidly and maximally rotating initial 
models Rl and R2, and the differentially rotating models 
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. — EOS with e* pairs 
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0.005 




0.02 
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t/M 




630 640 650 660 670 680 690 700 710 720 730 
t/M 

Fig. 2. — The upper panel displays the time evolution of the 
central rest-mass density for model Dl (a differentially rotating star 
with a mass M = 5 X 10 5 M Q and zero metallicity) with a T-law 
and the microphysical EOS with electron-positron pair creation. 
The middle and lower panels display the AH mass and the disk 
mass as a function of time for the collapse simulation with a T-law 
EOS. 

Dl and D2 are cconstructed with a polytropic EOS with 
the Lorene code (URL http:/ /www.lorene.obspm.fr). We 
obtain temperatures for our microphyscial models by in- 
verting the corresponding energy density with our EOS 
of Section 3.2. We also introduce a perturbation to trig- 
ger the gravitational collapse by reducing the pressure 
overall by « 1.5%. 

In order to determine the threshold metallicity re- 
quired to halt the collapse and produce an explosion we 
carry out several numerical simulations for each initial 
model with different values of the initial metallicity. The 
initial metallicities along with the fate of the star are 
given in Table 1. 

6. COMPARISONS WITH PREVtOUS STUDIES 
6.1. Comparison with ID calculations 

Axisymmctric calculations without rotation (i.e. mod- 
els SI and S2) retain the spherical symmetry of the initial 
conditions. There are no physical phenomena like con- 
vective or overturn instabilities to produce asphericity, 
i.e. we can directly compare our 2D non-rotating models 
with those computed in spher ical s ymmetry (ID calc u- 
lations) by IFuller et all (| 19861) and lLinke etaLI (|200lD . 

The mai n differences with r espect to the results ob- 
tained by IFuller et al.l (|1986[ ) are most likely due to 
two reasons. First, we apply a fully general relativis- 
tic treatment while they used a post-Newtonian treat- 
ment of gravity, and second, there are differences in the 
treatment of nuclear burning (Fuller et al. 1986 solved 
the relevant nuclear network without the approximations 
adopted in our work, see Section 3.3 for further details). 
Despite these differences, the results agree fairly well. 

2 In a core-collapse supernova nonradial instabilities are trig- 
gered either by negative entropy gradients caused by the shock 
deceleration and neutrino heating or by a generic instability of the 
stalled shock (SASI, Blondin and Mezzacapa 2003). Conditions for 
both processes are absent in the collapse and explosion of SMS. 
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Fig. 3. — Time evolution of the central rest-mass density for 
model D2 for three EOS, which shows that D2 becomes unsta- 
ble and collapses to a BH only when the microphysical EOS with 
electron-positron pairs is used. 

As discussed in detail in Section 7.1 the initial metallic- 
ities required to produce an explosion are similar, and a 
thermal bounce can be produced only if sufficient energy 
is liberated during the phase when the HCNO cycle is 
active. 

The main diffe rence with respect to the work of 
iLinke et al.l (|2001l) resides in the formulation of Ein- 
stein's field equations; in particular, in the foliation of 
the spacetime (foliation into a set spacelike hypersurfaces 
versus a foliation with outgoing null hypersurfaces) . I n 
order to compare with the results of ILinke et al.l |2001), 
we computed the redshifted total energy output for a 
model having the same rest-mass, doing the time inte- 
grati on until app r oxima tely the same evolutionary stage 
as in ILinke et al.l ([200 1). We find that the total energies 
released in neutrinos differ by less than 10% (for mode 
details see Section 7.3). 

6.2. T-law vs. microphysical EOS in uniformly rotating 

SMS 

Previous simulations of SMS collapse to BH in general 
relativity have been performed with a L-law EOS with 
r = 4/3 (with the only exception being the work of Linkc 
ct al. 2001). In order to elucidate the influence of the 
EOS on the dynamics of collapsing SMS, we performed 
three simulations of the same initial model (model R1.0, 
a marginally unstable uniformly rotating SMS with zero 
initial metallicity) without nuclear burning effects, and 
with three EOS: a T-law EOS with V = 4/3 (i.e. a similar 
set-up as in Shibata & Shapiro 2002) and the microphys- 
ical EOS with and without including electrons and the 
e pairs, i.e. for the last EOS case we consider Eqs. (19) 
and (20) for the baryons plus e 7 = aT 4 and P 1 — -|e 7 for 
the photons (where a is the radiation density constant). 
We denote the T-EOS as EOS-0, the full microphysical 
EOS, our canonical one for the studies of this work as 
EOS-1, and the reduced microphysical case as EOS-2. 

In Figure Q] we show with a dotted line the time evo- 
lution of the central density of model R1.0 with EOS-0, 
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Fig. 4. — Left upper panel shows the time evolution of the central rest-mass density for models SI and Rl (i.e., spherical and rotating 
stars with mass M = 5 X 10 5 Mq), and the lower left panel shows the time evolution of the central temperature. Horizontal dotted lines 
mark the temperature range in which nuclear energy is primarily released by the hot CNO cycle. Similarly, the time evolution of the 
same quantities for model S2 and R2 (i.e., spherical and rotating stars with mass M = 1 X 10 6 Mq) are shown in the upper and lower 
right panels. As the collapse proceeds, the central density and temperature rise rapidly, increasing the nuclear energy generation rate by 
hydrogen burning. If the metallicity is sufficiently high, enough energy can be liberated to produce a thermal bounce. This is the case for 
models Sl.c, Rl.d, S2.b and R2.d shown here. 



and with a dashed (solid) line the time evolution of the 
central density with EOS-2 (EOS-1). The first thing to 
note is that the collapse timescale obtained with the T- 
law EOS is shorter than that obtained with the micro- 
physical EOS (both EOS-1 and EOS-2), because the ion 
pressure contribution to the EOS raises the adiabatic in- 
dex above 4/3 (see Eq. [TB|) . This increase in the adiabatic 
index helps to stabilize the star against the gravitational 
instability, and therefore delays the collapse. 

On the other hand, the effect of pair creation reduces 
the adiabatic index below 4/3 at T > 10 9 K. This ex- 
plains the differences between the solid and dashed lines 
in Figure Q] at central densities p c > 10 g/cm 3 , which 
correspond to central temperatures T c > 10 9 K. Once the 
collapse enters this regime, pair creation becomes rele- 
vant enough to reduce the adiabatic index below 4/3, 
which destabilizes the collapsing star. Compared to pre- 
vious works (Shibata & Shapiro 2002) the use of a micro- 
physical EOS instead of a T-law EOS delays the collapse 
(mostly due to the baryons while destabilize) of an 
initially gravitationally unstable configuration. 

6.3. T-law vs. microphysical EOS in differentially 
rotating SMS 

We also performed 2D axisymmetric simulations of dif- 
ferentially rotating SMS. First, we investigated the influ- 
ence of the EOS on gravitationally unstable stars using 
model Dl as a reference. This model corresponds, within 
the accuracy to which the initial conditions can be repro- 
duc ed, to a different i ally ro tating unstable SMS discussed 
by iSaijo fe Hawkel (|2009f ) (i.e. their model I). Results 
are displayed in Figure [2] In the upper panel of this fig- 
ure, we show the time evolution of the central density for 
model Dl with EOS-0 (dashed line), and with the micro- 
physical EOS-1 (solid line). Opposite to the behavior in 
the case of the uniformly rotating SMS R1.0, the collapse 
timescale is longer with the T-law than with the micro- 



physical EOS. The reason for this difference is that the 
initial central temperature ( T c « 1.4 x 10 9 K) of model 
Dl is an order of magnitude higher than the initial cen- 
tral temperature in R1.0. Therefore, electron-positron 
pair creation reduces the stability of the star (by reducing 
T) already during the initial stages of the collapse. This 
behavior is expected to be present also in 3D, and since 
the collapse timescale is reduced when using the micro- 
physical EOS, nonaxisymmetric instabilities would have 
even less time to grow befor e the formation of a B H. It 
reinforces the conclusions of Saijo & Hawkc (2009), who 
showed that the three dimensional collapse of rotating 
stars proceeds in an approximately axisymmetric man- 
ner. 

The lower two panels of Figure [2] display the growth of 
the AH mass and the disk mass (defined as the rest-mass 
outside the AH of the newly formed BH) as a function 
of time (in units of the gravitational mass for compari- 
son with Figures 9 and 10 of Saijo & Hawke 2009), re- 
spectively. The values of both quantities at the end of 
the simulation agr ee, within a 5% differe nce, with those 
obtained in 3D bv ISaiio fc Hawkel (|2009f l. We note that 
there is also good agreement (less than 5% difference) re- 
garding the time at which an AH is first detected. These 
observed small differences are likely due to differences in 
the initial models and numerical techniques rather than 
to the influence of nonaxisymmtric effects. This suggests 
that our collapse simulation with the same treatment of 
ph ysics yields good agre ement with the 3D simulations 
of ISaiio fc Hawkel ([2009h 

We also investigated the influence of electron-positron 
pair creation on the evolution of gravitationally stable 
differentially rotating SMS using model D2 which is sim- 
ilar to the stable di fferentially rotating model III of 
ISaijo fc Hawkel (|2009l ). We performed three simulations 
of model D2 varying the EOS. In Figure [3] we show 
the central rest-mass density as a function of time for 
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a T-law (dashed line), and for the microphysical EOS- 
1 (solid line) and EOS-2 (d otted line). In agreem ent 
with the results obtained by iSaiio fc Hawkd (|2009f) we 
find that model D2 represents a stable differentially ro- 
tating SMS when a T-law is used. A persistent series 
of oscillations is triggered by the initial perturbation in 
the pressure. This is also the case with the microphysi- 
cal EOS-2 without the inclusion of electrons and the e ± 
pairs. However, the time evolution of D2 is completely 
different when pairs are taken into account. The in- 
fluence of pairs is large enough to destabilize model D2 
against gravitational collapse. We note that unlike all 
other SMS considered in this paper, which are T = 4/3 
models initially unstable to gravitational collapse, model 
D2 is an initially stable r = 4/3 model which becomes 
gravitationally unstable only by the creation of electron- 
positron pairs at high temperatures. Hence, using a mi- 
crophysical EOS with electron-positron pairs is crucial 
to determine the stability of differentially rotating SMS. 

We note that the central temperature of the initial 
models Dl and D2 is of the order of « 10 9 K. At this 
temperature, the main source of thermonuclear energy 
is hydrogen burning via the rp-process. It is however 
expected that such SMS would previously experience a 
phase of hydrogen burning via the cold and hot CNO 
cycles which would significantly affect the evolution of 
the models such that configurations with high T c as in 
models Dl and D2 might never be reached. Therefore, 
models Dl and D2 are not particularly well suited to 
investigate the existence of a thermal bounce during col- 
lapse (see Section 7). Exploring in detail the parameter 
space for the stability of differentially rotating SMS with 
the microphysical EOS, and the existence of a thermal 
bounce during the collapse phase depending on the ini- 
tial stellar metallicity, is a major task on its own, which 
is beyond the scope of this paper. For these reasons we 
do not consider differentially rotating SMS in this work. 

7. RESULTS 

7.1. Collapse to BH vs. Thermonuclear explosion 

First we consider a gravitationally unstable spheri- 
cally symmetric SMS with a gravitational mass of M = 
5 x lCrM© (SI. a, Sl.b and Six), which corresponds t o 
a model extensively discussed in iFuller et al.l (p.986), 
and therefore a llows for a compar ison with the results 
presented here. IFuller et al.l (|1986l ) found that unstable 
spherical SMS with M — 5 x 1O 5 M and an initial metal- 
licity Zcno = 2 x 10~ 3 collapse to a BH while models 
with an initial metallicity Zcno — 5 x 10 -3 explode due 
to the nuclear energy released by the hot CNO burning. 
They also found that the central density and tempera- 
ture at thermal bounce (where the collapse is reversed to 
an explosion) are p c ^ = 3.16 g/cm 3 and T c ^ = 2.6 x 10 8 
K, respectively. 

The left panels in Figure @] show the time evolution of 
the central rest-mass density (upper panel) and central 
temperature (lower panel) for models SI. a, Sl.c, Rl.a 
and Rl.d, i.e., non- rotating and rotating models with a 
mass of M = 5 x 10 5 Mq. In particular, the solid lines 
represent the time evolution of the central density and 
temperature for model Sl.c (Zcno — 7x 10~ 3 ) and Rl.d 
(Zcno = 5x 10~ 4 ). As the collapse proceeds, the cen- 
tral density and temperature rise rapidly, which increases 



the nuclear energy generation rate by hydrogen burning. 
Since the metallicity is sufficiently high, enough energy 
can be liberated to increase the pressure and to produce a 
thermal bounce. This is the case for model Sl.c. In Fig- 
ure |4] we show that a thermal bounce occurs (at approxi- 
mately t ~ 7 x 10 5 s) entirely due to the hot CNO cycle, 
which is the main source of thermonuclear energy at tem- 
peratures in the range 2 x 10 8 K < T < 5 x 10 8 K. The 
rest-mass density at bounce is p c .\, — 4.8 g/cm 3 and the 
temperature T c ,b = 3.05 x 10 8 K. These values, as well as 
the threshold metallicity needed to trigger a thermonu- 
clear expl osion (Zcno = 7x 1 0~ 3 ), are higher than those 
found by IFuller et al.l (|1986l ) (who found that a spher- 
ical nonrotating model with the same rest-mass would 
explode, if the initial metallicity was Zcno = 5x 10 -3 ). 

On the other hand, dashed lines show the time evo- 
lution of the central density and temperature for model 
SI. a (Zcno = 5 x 10~ 3 ). In this case, as well as for 
model Sl.b, the collapse is not halted by the energy re- 
lease and continues until an AH is found, indicating the 
formation of a BH. 

We note that the radial velocity profiles change contin- 
uously near the time where the collapse is reversed to an 
explosion due to the nuclear energy released by the hot 
CNO burning, and an expanding shock forms only near 
the surface of the star at a radius R « 1.365 x 10 13 cm 
(i.e. R/M w 180) where the rest-mass density is « 
3.5 x 10 _6 gcm -3 . We show in Figure[5]the profiles of the 
^-component of the three-velocity v x along the x-axis 
(in the equatorial plane) for the nonrotating spherical 
stars SI. a (dashed lines) and Sl.c (solid lines) at three 
different time slices near the time at which model Sl.c 
experiences a thermal bounce. Velocity profiles of model 

51. c are displayed up to the radius where a shock forms 
at t w 7.31 x 10 5 s and begins to expand into the low 
density outer layers of the SMS. 

The evolutionary tracks for the central density and 
temperature of the rotating models Rl.a and Rl.d are 
also shown in Figure [4] A dashed line corresponds to 
model Rl.a, with an initial metallicity Zcno = 5 x 10~ 4 , 
which collapses to a BH. A solid line denotes model Rl.d 
with Zcno — 2 x 10 -3 , which explodes due to the energy 
released by the hot CNO cycle. We find that Model Rl.c 
with a lower metallicity of Zc no = 1 X 10 -3 also explodes 
when the central temperature is the range dominated by 
the hot CNO cycle. 

As a result of the kinetic energy stored in the rotation 
of models Rl.c and Rl.d, the critical metallicity needed 
to trigger an explosion decreases significantly relative to 
the non-rotating case. We observe that rotating models 
with initial metallicities up to Zcno = 8 x 10~ 4 do not 
explode even via the rp-process, which is dominant at 
temperatures above T m 5 x 10 8 K and increases the 
hydrogen burning rate by 200 — 300 times relative to 
the hot CNO cycle. We also note that the evolution 
time scales of the collapse and bounce phases are reduced 
because rotating models are more compact and have a 
higher initial central density and temperature than the 
spherical ones at the onset of the gravitational instability. 

The right panels in Figure [4] show the time evolu- 
tion of the central rest-mass density (upper panel) and 
the central temperature (lower panel) for models S2.a, 

52. b, R2.a and R2.d, i.e., of models with a mass of 
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Fig. 5. — Profiles of the x-componcnt of the three-velocity v x 
along the ir-axis (in the equatorial plane) for the nonrotating spher- 
ical stars SI. a (dashed lines) and Sl.c (solid lines) at three different 
time slices near the time at which model Sl.c experiences a ther- 
mal bounce. Velocity profiles of model Sl.c are displayed up to 
the radius where a shock, that expands into the low density outer 
layers of the SMS, forms. 



M = 10 6 M©. We find that the critical metallicity for an 
explosion in the spherical case is Zqno = 5x 10 -2 (model 
S2.b), while model S2.a with Zcno = 3 x 10~ 2 collapses 
to a BH. We note that the critical metallicity leading to a 
thermon uclear explos i on is h igher than the critical value 
found bv lFuller et"aH (|1986l) (i.e., Z CNO = 1 X 10~ 2 ) for 
a spherical SMS with the same mass. The initial metal- 
licity leading to an explosion in the rotating case (model 
R2.d) is more than an order of magnitude smaller than in 
the spherical case. As for the models with a smaller grav- 
itational mass, the thermal bounce takes place when the 
physical conditions in the central region of the star allow 
for the release of energy by hydrogen burning through 
the hot CNO cycle. Overall, the dynamics of the more 
massive models indicates that the critical initial metal- 
licity required to produce an explosion increases with the 
rest-mass of the star. 

Figure [5] shows the total nuclear energy generation rate 
in erg/s for the exploding models as a function of time 
during the late stages of the collapse just before and af- 
ter bounce. The main contribution to the nuclear energy 
generation is due to hydrogen burning by the hot CNO 
cycle. The peak values of the energy generation rate 
at bounce lie between several 10 51 erg/s for the rotating 
models (Rl.d and R2.d), and « 10 52 - 10 53 erg/s for the 
spherical models (Sl.c and S2.b). As expected the max- 
imum nuclear energy generation rate needed to produce 
an explosion is lower in the rotating models. Moreover, 
as the explosions are due to the energy release by hy- 
drogen burning via the hot CNO cycle, the ejecta would 
mostly be composed of 4 He. 

As a result of the thermal bounce, the kinetic energy 
rises until most of the energy of the explosion is in the 
form of kinetic energy. We list in the second but last 
column of Table 1 the radial kinetic energy after thermal 
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Fig. 6. — Nuclear energy generation rate in erg/s for the explod- 
ing models (Sl.c, Rl.d, Rl.c, S2.b and R2.d) as a function of time 
near the bounce. The contribution to the nuclear energy gener- 
ation is mainly due to hydrogen burning by the hot CNO cycle. 
The peak values of the energy generation rate at bounce lie be- 
tween 10 [erg/s] for the rotating models (Rl.d and R2.d), and 



10° 



10 [erg/s] for the spherical models (Sl.c and S2.b). 



bounce, which ranges between -Erk = 1-0 x 10 55 ergs for 
the rotating star Rl.c, and £"rk = 3.5 x 10 57 ergs for the 
spherical star S2.b. 

7.2. Photon luminosity 

Due to the lack of resolution at the surface of the star, 
it becomes difficult to compute accurately the photo- 
sphere and its effective temperature from the criterion 
that the optical depth is r = 2/3. Therefore, in order to 
estimate the photon luminosity produced in association 
with the thermonuclear explosion, we make use of the 
fact that within the diffusion approximation the radia- 
tion flux is given by 



3n es p 



VU, 



(35) 



where U is the energy density of the radiation, and n es is 
the opacity due to electron Thompson scattering, which 
is the main source of opacity in SMS. The photon lumi- 
nosity in terms of the temperature gradient and for the 
spherically symmetric case can be written as 



167racr 2 T 3 dT 
3n es p dr ' 



(36) 



where a is the radiation constant, and c the speed of light. 
As can be seen in the last panel of Figure [7] the distribu- 
tion of matter becomes spherically symmetric during the 
phase of expansion after the thermonuclear explosion. In 
this figure (Fig. [7]) we show the isodensity contours for 
the rotating model Rl.d. The frames have been taken at 
the initial time (left figure), at t = 0.83 x 10 5 s (central 
figure) just after the thermal bounce (at t = 0.78 x 10 5 s), 
and at t = 2.0 x 10 5 s when the radius of the expanding 
matter is roughly 4 times the radius of the star at the 
onset of the collapse. 
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Fig. 7. — Isodensity contours of the logarithm of the rest-mass density (in g/cm 3 ) for the rotating model Rl.d. The frames have been 
taken at the initial time (left figure), at t = 0.83 X 10 5 s (central figure) just after a thermal bounce takes place, and at t = 2.0 X 10 s 
when the radius of the expanding matter is roughly 4 times the radius of the star at the onset of the collapse. 

ity iFuller et all (|1986l) found for a nonrotating SMS of 
same rest-mass. The photon luminosity remains super- 
Eddington during the phase of rapid expansion that fol- 
lows the thermal bounce. We compute the photon lu- 
minosity until the surface of the star reaches the outer 
boundary of the computational domain « 1.0 x 10 5 
after the bounce. Beyond that point, the luminosity 
is expected to decrease, and then rise to a plateau of 
~ 10 45 erg/s due to the recombination of hydrogen (see 
Fuller et al. 1986 for a nonrotating star). 

7.3. Collapse to BH and neutrino emission 

The outcome of the evolution of models that do not 
generate enough nuclear energy during the contraction 
phase to halt the collapse is the formation of a BH. 
The evolutionary tracks for the central density and tem- 
perature of some of these models are also shown in 
Figure 0J The central density typically increases up 
to p c ~ 10 7 gcm~ 3 and the central temperature up to 
T c ~ 10 10 K just before the formation of an AH. 

Three isodensity contours for the rotating model Rl.a 
collapsing to a BH are shown in Figure [9l which display 
the flattening of the star as the collapse proceeds. The 
frames have been taken at the initial time (left panel), 
at t = 0.83 x 10 5 s (central panel) approximately when 
model Rl.d with higher metallicity experiences a thermal 
bounce and, at t = 1.127 x 10 5 s, where a BH has already 
formed and its AH has a mass of 50% of the total initial 
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Fig. 8. — Logarithm of the photon luminosity of model Rl.d 
in units of erg/s as a function of time. The vertical dashed line 
indicates the time at which the thermal bounce takes place. 

The photon luminosity computed using Eq. (|36[) for 
model Rl.d is displayed in Figure [SJ where we also in- 
dicate with a dashed vertical line the time at which the 
thermal bounce takes place. The photon luminosity be- 
fore the thermal bounce is computed at radii inside the 
star unaffected by the local dynamics of the low den- 
sity outer layers which is caused by the initial pressure 
perturbation and by the interaction between the surface 
of the SMS and the artificial atmosphere. Once the ex- 
panding shock forms near the surface, the photon lumi- 
nosity is computed near the surface of the star. The 
lightcurve shows that, during the initial phase, the lu- 
minosity is roughly equal to the Eddington luminosity 
« 5 x 10 43 erg/s until the thermal bounce. Then, the 
photon luminosity becomes super-Eddington when the 
expanding shock reaches the outer layers of the star and 
reaches a value of L 1 « 1 x 10 45 erg/s. This value of 
the photon luminosity after the bounce is within a few 
percent difference with respect to the photon luminos- 



At the temperatures reached during the late stages of 
the gravitational collapse (in fact at T > 5 x 10 8 K) the 
most efficient process for hydrogen burning is the break- 
out from the hot CNO cycle via the 15 0(a, 7) 19 Ne re- 
action. Nevertheless, we find that models which do not 
release enough nuclear energy by the hot CNO cycle to 
halt their collapse to a BH, are not able to produce a 
thermal explosion due to the energy liberated by the 
15 0(a,7) 19 Ne reaction. We note that above 10 9 K, not 
all the liberated energy is used to increase the temper- 
ature and pressure, but is partially used to create the 
rest-mass of the electron-positron pairs. As a result of 
pair creation, the adiabatic index of the star decreases, 
which means the stability of the star is reduced. More- 
over, due to the presence of e ± pairs, neutrino energy 
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Fig. 9. — Isodensity contours of the logarithm of the rest-mass density (in g/cm 3 ) for the rotating model Rl.a. The frames have been 
taken at the initial time (left panel), at t = 0.83 X 10 5 s (central panel), at t = 1.127 X 10 5 s, where a BH has already formed and its 
apparent horizon encloses a mass of 50% of the total initial gravitational mass. 

losses grow dramatically. 

Figure [10] shows (solid lines) the time evolution of the 
rcdshifted neutrino luminosities of four models collaps- 
ing to a BH (SI. a, Rl.a, S2.a, and R2.a), and (dashed 
lines) of four models experiencing a thermal bounce 
(Sl.c, Rl.d, S2.b, and R2.d). The change of the slope 
of the neutrino luminosities at ~ 10 43 erg/s denotes the 
transition from photo-neutrino emission to the pair an- 
nihilation dominated region. The peak luminosities in 
all form of neutrino for models collapsing to a BH are 
L v ~ 10 55 erg/s. Neutrino luminosities can be that im- 
portant because the densities in the core prior to BH 
formation are p c ~ 10 7 gcm~ 3 , and therefore neutrinos 
can escape. The peak ne utrino luminosities lie between 
the luminosities found by iLinke et al.l (|2001|) for the col- 
lapse o f spherical SMS, and those found bv lWooslev et al.l 
( 1986) (who only took into account the luminosity in the 
form of electron antineutrino) . The maximum luminos- 
ity decreases slightly as the rest-mass of the i nitial model 
increa ses, which was already observed by ILinke et al.1 
(|2001h . In addition, we find that the peak of the red- 
shifted neutrino luminosity does not seem to be very sen- 
sitive to the initial rotation rate of the star. We also note 
that the luminosity of model Rl.a reflects the effects of 
hydrogen burning at L v ~ 10 43 erg/s. 

The total energy output in the form of neutrinos is 
listed in the last column of Table 1 for several models. 
The total radiated energies vary between E v ~ 10 56 ergs 
for models collapsing to a BH, and E v ~ 10 45 — 10 46 
ergs for exploding models. These results are in reason- 
able agreement with previous calculations. For instance, 
iWooslev et al.l (|1986l ) obtained that the total energy out- 
put in the form of electron antineutrinos for a spherical 
SMS with a mass 5 x 10 5 M Q and zero initial metallic- 
ity was 2.6 x 10 56 ergs, although their simulations ne- 
glected general relativistic effects which are important 
to compute accurately the rela tivistic redshifts. On the 
other hand. ILinke et all (j2001l ). by means of relativistic 
one-dimensional simulations, found a total radiated en- 
ergy in form of neutrinos of about 3 x 10 56 ergs for the 
same initial model, and about 1 x 10 56 ergs when red- 
shifts were ta ken into account. I n order to compare with 
the results of ILinke et a l. (2001), we computed the red- 
shifted total energy output for Model SI. a, having the 




Fig. 10. — Time evolution of the redshifted neutrino luminosities 
for models Rl.a, SI. a, R2.a, and S2.a all collapsing to a BH; and for 
models Rl.d, Sl.c, R2.d, and S2.b experiencing a thermal bounce. 
The time is measured relative to the collapse timescale of each 
model, Rl, SI, R2 and S2, with t (1, 7, 0.2, 6) in units of 10 5 s. 

same res t-mass, unti l appro ximately the same evolution 
stage as ILinke et al.l (|2001[) did (i.e. when the differen- 
tial neutrino luminosity dL v /dr ~ 4 x 10 45 erg/s/cm). 
We find that the total energy released in neutrinos is 
1.1 x 10 56 ergs. 

The neutrino luminosities for models experiencing a 
thermonuclear explosion (dashed lines in Fig. ITU)) peak 
at much lower values L v ~ 10 42 — 10 43 erg/s, and decrease 
due to the expansion and disruption of the star after the 
bounce. 

7.4. Implications for gravitational wave emission 

The axisymmetric gravitational collapse of rotating 
SMS with uniform rot ation is expected to emit a burst 
of gr avitational waves (jSaiio et al.l(200l ISaiio fc Hawkel 
2009) with a frequency within the LISA low frequency 
band (10~ 4 — 10 _1 Hz). Although through the Simula- 
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tions presented here we could not investigate the devel- 
opment of nonaxisymmetric features in our axisymmet- 
ric models that could also lead to the emission of GWs, 
Saiio & Hawkej (|2009P) have shown that the three dimen- 
sional collapse of rotating stars proceeds in an approxi- 
mately axisymmetric manner. 

In an axisymmetric spacetime, the x-mode vanishes 
and the +- mode of gravitational waves with I = 2 
computed using the quadr upole formula is written as 
(jShibata fc Sek iguchi 2003) 



io- 19 



i quad Ixx{tret) ^-zzi^ret) 

= SIT 



(37) 



where Iij refers to the second time derivative of the 
quadrupole moment. The gravitational wave quadrupole 
amplitude is A 2 (t) = I xx (t ret ) - I Z z(Uet)- Following 
IShibata fc Sekiguchl (|2003l ) we compute the second time 
derivative of the quadrupole moment by finite differenc- 
ing the numerical results for the first time derivative of 
hj obtained by 



hj = I P* (v l x j + x t v j ) d 3 x. 



(38) 



We calculate the characteris tic gravitational wave 
strain (Flanagan & Hughes 1998) as 



•(/) 



12 G 1 dE{f) 
n~c?rP df : 



(39) 



where D is the distance of the source, and dE(f)/df 
the spectral energy density of the gravitational radiation 
given by 

dE(f) _ c 3 (2^/) 2 - 2 



with 



df G 167T Mf) 
Mf) = / A 2 {t)e 2 ^ l dt 



(40) 
(41) 



We have calculated the quadrupole gravitational wave 
emission for the rotating model Rl.a collapsing to a BH. 
We plot in Figurc[TT]thc characteristic gravitational wave 
strain (Eq |39[) for this model assuming that the source 
is located at a distance of 50 Gpc (i.e., z » 11) , to- 
gether with the desig n noise spectrum h (f ) = yfStJj) 



of the LISA detecto r (lLarson et al.ll2000D We find that 
in ag reeme nt with ISaiio et al. ( 20021 ) . ISaiio fc Hawkd 
(2009) and iFrver fc New! (|2011f ). the burst of gravita- 
tional waves due to the collapse of a rotating SMS could 
be detected at a distance of 50 Gpc and at a frequen cy 
which approximately takes the form (|Saiio et alj r2002) 



fi 



burst 



3 x 10" 



10 b M, 



/5M 



3/2 



[Hzl 



(42) 



where R/M is a characteristic mean radius during black 
hole formation ( typically set t o R/M = 5). 

Furthermore, IKiuchi et al.l (|2011[) have recently in- 
vestigated, by means of three-dimensional general rel- 
ativistic numerical simulations of equilibrium tori or- 
biting BHs, the development of the nonaxisymmetric 
Papaloizou-Pringle instability (PPI) in such systems 




o.oooi 



Fig. 11. — Characteristic gravitational wave strain for model Rl.a 
assuming that the source is located at a distance of 50 Gpc, to- 
gether with the design noise spectrum h(f) = <J fSfo (/) for LISA 
detector. 



(Papaloizou & Pringlc 1984), and have found that a non- 
axisymmetric instability associated with the m — 1 mode 
grows for a wide range of self-gravitating tori orbiting 
BHs, leadi ng to the emission of quasiperiodic GWs. In 
particular, IKiuchi et al.l ()2011[ ) have pointed out that the 
emission of quasiperiodic GWs from the torus resulting 
after the formation of a SMBH via the collapse of a SMS 
could be well above the noise sensitivity curve of LISA for 
sources located at a distance of lOGpc. Such instability 
appears for tori whose angular velocity in the equato- 
rial plane expressed as H(r) oc r q has q < q^ep where 
qkep corresponds to the Keplerian limit, i.e. q — — 1.5 in 
Newtonian gravity. 

We find that the torus (defined as the rest-mass out- 
side the AH) that forms after the collapse to a BH of the 
uniformly rotating model Rl.a (when the mass of the AH 
exceeds 50% of the ADM mass) does not fulfill the above 
condition for the development of the PPI. However, we 
find that the torus that forms when the differentially ro- 
tating model Dl collapses to a_BH has a distribution of 
angular momentum such that Q(r) oc r q with q « —1.62. 
This suggests that the torus may be prone to the devel- 
opment of the nonaxisymmetric PPI, which would lead 
to the emission of quasiperiodic GWs with peak ampli- 
tude ~ 10~ 18 — 10~ 19 and frequency ~ 10 _3 Hz during 
an accretion timescale ~ 10 5 s. 

7.5. Conclusions 

We have presented results of general relativistic simu- 
lations of collapsing supermassive stars using the two- 
dimensional general relativistic numerical code Nada, 
which solves the Einstein equations written in the BSSN 
formalism and the general relativistic hydrodynamic 
equations with high resolution shock capturing schemes. 
These numerical simulations have used an EOS that in- 
cludes the effects of gas pressure, and tabulated those 
associated with radiation pressure and electron-positron 
pairs. We have also taken into account the effects of 
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thermonuclear energy release by hydrogen and helium 
burning. In particular, we have investigated the effects 
of hydrogen burning by the ^-limited hot CNO cycle and 
its breakout via the 15 0(a, 7) 19 Ne reaction (rp-process) 
on the gravitational collapse of nonrotating and rotating 
SMS with non-zero metallicity. 

We have presented a comparison with previous studies, 
and investigated the influence of the EOS on the collapse. 
We emphasize that axisymmetric calculations without 
rotation (i.e. models SI and S2) retain the spherical sym- 
metry of the initial configurations as there are no physi- 
cal phenomena to produce asphericity and numerical ar- 
tifacts associated with the use of Cartesian coordinates 
are negligibly small. Overall, our collapse simulations 
yield good agreement with previous works when using 
the same treatment of physics. We have also found that 
the collapse timescale depends on the ion contributions 
to the EOS, and electron-positron pair creation affects 
the stability of SMS. Interestingly, differentially rotating 
stars that are gravitationally stable with a V = 4/3 EOS 
can become unstable against gravitational collapse when 
the calculation is performed with the microphysical EOS 
including pair creation. 

We have found that objects with a mass of « 5 x 10 5 M Q 
and an initial metallicity greater than Zqno ~ 0.007 
explode if non-rotating, while the threshold metallicity 
for an explosion is reduced to Zcno ~ 0.001 for ob- 
jects which are uniformly rotating. The critical initial 
metallicity for a thermal explosion increases for stars 
with a mass of w 10 6 M Q . The most important con- 
tribution to the nuclear energy generation is due to the 
hot CNO cycle. The peak values of the nuclear energy 
generation rate at bounce range from ~ 10 51 erg/s for ro- 
tating models (Rl.d and R2.d), to - 10 52 - 10 53 erg/s 
for spherical models (Six and S2.b). After the ther- 
mal bounce, the radial kinetic energy of the explosion 



rises until most of the energy is kinetic, with values 
ranging from Ek ~ 10 56 ergs for rotating stars, to up 
Ek ~ 10 57 ergs for the spherical star S2.b. The neutrino 
luminosities for models experiencing a thermal bounce 
peak at L v ~ 10 42 erg/s. 

The photon luminosity roughly equal to the Edding- 
ton luminosity during the initial phase of contraction. 
Then, after the thermal bounce, the photon luminosity 
becomes super-Eddington with a value of about L 1 « 
1 x 10 45 erg/s during the phase of rapid expansion that fol- 
lows the thermal bounce. For those stars that do not ex- 
plode we have followed the evolution beyond the phase of 
black hole formation and computed the neutrino energy 
loss. The peak neutrino luminosities are L v ~ 10 55 erg/s. 

SMS with masses less than w 10 6 M Q could have 
formed in massive halos with T„j r > 10 4 K. Although the 
amount of metals that was present in such environments 
at the time when SMS might have formed is unclear, 
it seems possible that the metallicities could have been 
smaller than the critical metallicities required to reverse 
the gravitational collapse of a SMS into an explosion. If 
so, the final fate of the gravitational collapse of rotating 
SMS would be the formation of a SMBH and a torus. In 
a follow-up paper, we aim to investigate in detail the dy- 
namics of such systems (collapsing of SMS to a BH-torus 
system) in 3D, focusing on the post-BH evolution and 
the development of nonaxisymmetric features that could 
emit detectable gravitational radiation. 
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Forschungsgesellschaft (DFG) through its Transregional 
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